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£^ I Abstract. We present a general and practical procedure to solve the general relativistic 
hydrodynamic equations by using any of the special relativistic Riemann solvers recently 
OO . developed for describing the evolution of special relativistic flows. Our proposal relies on a 
local change of coordinates in terms of which the spacetime metric is locally Minkowskian 
Oh^ and permits accurate numerical calculations of general relativistic hydrodynamics prob- 
Q [ lems using the numerical tools developed for the special relativistic case with negligible 

r/3 • computational cost. The feasibility of the method has been confirmed by a number of 

■ 

r* ■ numerical experiments. 
• i— i . 
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1. Introduction 

In the near future the first generation of Earth-based laser-interferometer detectors of 
gravitational waves will be operating (LIGO, VIRGO, GEO600, TAMA). This perspec- 
tive has stimulated researchers working in numerical relativistic astrophysics to develop 
robust codes for the simulation of the different astrophysical sources of gravitational ra- 
diation, such as, e.g., stellar core collapse, coalescing compact binaries or accretion onto 
compact objects. 

Relativistic hydrodynamical codes experienced a substantial advance at the beginning 



of the nineties (Marti, Ibahez, Miralles, 1991 ) with the implementation of high-resolution 
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shock-capturing methods (HRSC) originally developed in classical fluid dynamics. These 
methods make use of the hyperbolic and conservative character of the equations and dis- 
play a number of interesting features and properties, as being stable and conservative, 
converging to physical solutions and having high accuracy in regions where the solu- 
tion is smooth. They are all based on the resolution of local Riemann problems at the 



interfaces of numerical cells -following the seminal idea of Godunov (1959) - ensuring 
a consistent treatment of discontinuities (shocks). The first relativistic applications of 
these techniques showed their capabilities in describing accurately complex flows, with 



high Lorentz factors and strong shocks, superseding traditional methods (Wilson 197E) 



which failed to describe ultrarelativistic flows (Norman and Winkler 1986 ). Up to now, 
the use of HRSC methods in relativity has been mainly restricted to the field of spe- 
cial relativistic hydrodynamics (SRH) in the simulation of collisions of heavy ions and, 



remarkably, extragalactic jets (Marti, Muller and Ibanez, 1994). We refer the interested 



reader to the introductory section in Banyuls et al. (1997) for a recent review of the 
current status of HRSC techniques in numerical relativistic hydrodynamics. 

In recent years, the great success obtained in the first relativistic applications has 
drawn the attention of specialists who started to develop specific Riemann solvers for 
SRH. Nowadays, most of the reliable HRSC methods developed during the last twenty 
years in classical hydrodynamics have their special-relativistic extension (see, e.g., Ibanez 



et at, 1997, for an updated list). In the case of general relativistic hydrodynamics (GRH), 
the development of numerical codes based on the resolution of local Riemann problems is 
still in its infancy. Only a small number of papers have considered the extension of HRSC 
methods from SRH to GRH (see below). In addition, recently, several formulations of 
Einstein Equations as a first-order hyperbolic system of balance laws have been derived 



(Friedrich |1985| ; Bona et al. |1995| ; Abrahams et al. |1995| ; Fritelli & Reula |1996|). This 
opens a new strategy, in the field of Numerical Relativity, permitting the use of HRSC 
schemes, specifically designed for such hyperbolic systems, to solve the coupled system 



of spacetime plus hydrodynamics (Bona et al. 1993). 

The basic idea behind this work is to obtain a general procedure that allows us to 
take advantage of the increasing number of special relativistic Riemann solvers (SRRS) 
developed recently, in order to generate numerical solutions describing the evolution of 
relativistic flows in strong gravitational fields (background or dynamical) . All the previous 
works done to extend HRSC methods to GRH have used linearized Riemann solvers 



(Marti, Ibanez, Miralles 1991; Eulderink and Mellema 1994; Romero et al. 1996; Banyuls 



et al. 1997). In this paper we describe a procedure to use any type of SRRS in general 
relativistic hydrodynamics, including the exact Riemann solver for ID problems. This 
procedure relies on a local change of coordinates at each numerical interface, in terms of 
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which the spacetime metric is locally flat, analogously to the approach followed in classical 
fluid dynamics, when using the solution of Riemann problems in general curvilinear 
coordinates. The numerical implementation is simple, computationally inexpensive, and 
provides a useful tool for the researchers currently working in SRH to enter the field of 
GRH. 

The structure of this paper is the following: In §2 we summarize the GRH equations 
in the {3 + 1} formalism and the basic ideas of HRSC methods. In §3 the formulation 
of Riemann problems in locally flat spacetimes and the method to obtain the numer- 
ical fluxes for the GRH equations from those obtained in SRH is explained. In §4 we 
briefly describe the set of numerical tests and applications performed to demonstrate 
the feasibility of the approach. Finally, in §5 we summarize our results and foresee other 
applications of our proposal. 

2. GRH equations in the 3+1 formalism and HRSC methods. 

Let be a general spacetime, described by the four dimensional metric tensor g. The 
line element has the form 

ds 2 = -{a 2 - (3 t (3 l )dt 2 + 2(3 i dx i dt + j lj dx t dx j (1) 

Throughout the paper Greek (Latin) subscripts run from to 3 (1 to 3) and geometrized 
units are used (G = c = 1). Let {d t , di} be the coordinate basis and n the unit timelike 
vector field normal to the spacelike hypersurfaces E t (t = const.) 

dt = an + (3*d t (2) 

We denote by J and T the density current and the energy-momentum tensor for a 
perfect fluid, respectively, 

J = pu (3) 

T = phu <g> u + pg (4) 

with u the four-velocity of the fluid, p the rest-mass density, p the pressure and h the 
specific enthalpy (h = 1 + e + p/ p, where e is the specific internal energy). 

The equations describing the evolution of a relativistic fluid are the local conservation 
laws of baryon number and energy-momentum and can be written, for observers which 
are at rest in the slice S t (Eulerian observers), in terms of the divergence of the 5 vector 
fields {J, T • n, T • d u T • d 2 , T • d 3 } as, 

V-A = s, (5) 

where A denotes any of the above 5 vector fields and s is the corresponding source term. 
Explicit expressions of these vectors in terms of the primitive variables {p, e, v{\ (with Vi 
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the components of the velocity measured by an Eulerian observer) , as well as expressions 
for the source terms, are given in Banyuls et al. ( 1997] ). 



Let us consider the integral form of the above equations on a four-dimensional volume 
Q C M. with three-dimensional boundary c?il, and apply Gauss theorem to obtain the 
corresponding balance equation 

A-d 3 £ = [ sdQ. (6) 
an Jn 

For numerical applications, we choose the volume £1 as the one bounded by the co- 
ordinate hypersurfaces {S^q, ~S x a +Ax c "}- Hence, the time variation of the mean value of 
A , 

i rx 1 +Ax 1 r x 2 +Ax 2 r x 3 +Ax 3 

A° = _L / / / ^r g A»dxHx*dx\ (7) 



AV . . 

within the spatial volume 

px 1 +Ax 1 px 2 +Ax 2 px 3 +Ax 3 

AV= / / / ^^dx 1 dx 2 dx 3 , (8) 

J x 1 J x 2 J X 3, 

can be obtained from 



(A°AV) t+At = (A°AV)i + / sdQ 

Jn 



A-d 3 T,+ I A-d 3 £- 



A • d 3 E - / A • ,/ :; S 

3i 



A ■ d 23 + / A d A T, . (9) 



In order to update the solution in time, the volume and surface integrals on the right 
hand side of Eq (^|) have to be evaluated (see §3). HRSC schemes rely on the calculation 
of the A vector fields by solving local Riemann problems combined with monotonized 
cell reconstruction techniques. 

3. Formulation of Riemann Problems in locally Minkowskian coordinates. 

According to the equivalence principle, physical laws in a local inertial frame of a curved 
spacetime have the same form as in special relativity. Locally flat (or geodesic) systems of 
coordinates, in which the metric is brought into the Minkowskian form up to second order 
terms, are the realization of these local inertial frames. However, whereas the coordinate 
transformation leading to locally flat coordinates involves second order terms, locally 
Minkowskian coordinates are obtained by a linear transformation. This fact is of crucial 
importance when exploiting the selfsimilar character of the solution of the Riemann 
problems set up across the coordinate surfaces. 
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Hence, we propose to perform, at each numerical interface, a coordinate transfor- 
mation to locally Minkowskian coordinates assuming that the solution of the Riemann 
problem is the one in special relativity and planar symmetry. This last assumption is 
equivalent to the approach followed in classical fluid dynamics, when using the solution 
of Riemann problems in slab symmetry for problems in cylindrical or spherical coor- 
dinates, which breaks down near the singular points {e.g. the polar axis in cylindrical 
coordinates). Analogously to classical fluid dynamics, the accuracy will depend on the 
magnitude of the Christoffel symbols, which might be large whenever huge gradients or 
large temporal variations of the gravitational field are present. Finer grids and improved 
time advancing methods will be required in those regions. 

In the rest of this section we will focus on the evaluation of the first flux integral in 
Eq. (|). 

A-d 3 S= f A 1 x /^dx°dx 2 dx 3 (10) 

To begin, we will express the integral on a basis with eg = n and forming an 
orthonormal basis in the plane orthogonal to n with ej normal to the surface E s i and 
eg and tangent to that surface. The vectors of this basis verify • = r) & & with rj & 
the Minkowski metric (in the following, caret subscripts will refer to vector components 
in this basis). 

Denoting by Xq the coordinates of the center of the interface at time t, we introduce 
the following locally Minkowskian coordinate system 

x & = M«(x'*-x%), (11) 

where the matrix is given by d a = AZ^e^, calculated at Xq. In this system of 
coordinates the flux terms in the equations of GRH are written as in SRH, in Cartesian 
coordinates, and the flux integral ([h]) reads 

CA 1 - ^-A^J^jdx^dx^dx 3 (12) 
=.1 a 

with yf—g = 1 + 0(x a ), where we have taken into account that, in the coordinates x a , 
E a i is described by the equation x^ — @—x^ — (with f3 l = where the metric 

elements (3 1 and a are calculated at Xq. Therefore, the effect of a non-zero shift is that 
the interface where the Riemann problem has to be solved is not at rest but moves with 
speed 

At this point, all the theoretical work developed in the last years, concerning SRRS, 
can be exploited. The procedure involves the following steps: 

1) We set up the Riemann problem by giving the values at the two sides of Y, x i of 
two independent thermodynamical variables (namely, the rest mass density p and the 



6 Jose A. Pons et ai: GRH with SRRS 

specific internal energy e) and the components of the velocity in the orthonormal spatial 
basis v l , which are calculated using 



(13) 



where 



Ml 



-7 722 +7 723 

7 11 VT22 
"v/722 

723 



/722 



-7 13 V722733~(723) 
7 11 %/722 



V 722 733 -(723) 2 
V722 



(14) 



2) The special relativistic Riemman problem is solved for the variables p, e and v\ 
obtaining the fluxes associated to J, T • n, T • e- . Notice that the effect of a non-zero shift 
has to be considered at this stage. Although most linearized Riemann solvers provide 
the numerical fluxes for surfaces at rest, it is easy to generalize them to moving surfaces 
relying on the conservative and hyperbolic character of the system of equations (see, e. 



g., Harten and Hyman 1983 ). 

3) Once the Riemann problem has been solved, by means of any linearized or exact 
SRRS, we can take advantage of the selfsimilar character of the solution of the Rie- 
mann problem, which makes it constant on the surface E x i simplifying enormously the 
calculation of the above integral (|l2]): 



A • d 3 S = {A 1 - —A 8 ) 
a 



-g dx°dx 2 dx 3 



(15) 



where the superscript (*) stands for the value on T, x i obtained from the solution of the 
Riemann problem. The quantity in parenthesis in Eq. ( |l5| ) represents the numerical flux 
across S^i . Notice that the numerical fluxes correspond to the vector fields J, T • n, T • ej, 
T • eg, T • eg. In order to obtain the momentum fluxes in the original coordinates (T • di) 
the additional relation 



T.5 i = M/(T.e,) 



(16) 



has to be used. 

4) Finally, the numerical fluxes are multiplied by the surface integral appearing at 
the right hand side of (]lq), that is expressed in terms of the original coordinates as 



-g dx°dx 2 dx 3 



(17) 



and can be easily evaluated for a given metric. 



Let us remind that, in this section, we have focussed on the calculation of the flux 
terms in Eq. (^|), for given left and right states. Obviously, the performance of the numer- 
ical code depends on the quality of the provided initial states, as well as the computation 
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of the source terms in Eq. (^), and the algorithm for time advancing. In all the calcu- 
lations presented in next section, left and right states for Riemann problems have been 
obtained with a monotonic, piecewise linear reconstruction procedure. The source inte- 
grals have been evaluated assuming constant values of p, e and v l inside the numerical 
cells. Finally, advance in time has been done by means of a TVD-preserving, third order 
Runge-Kutta, that takes into account the influence of the source terms in the interme- 
diate steps. Notice that the treatment of the source terms is essential for the method 
to work in regions where they dominate. A treatment consistent with the reconstruction 
procedure and better time advancing schemes may be required in regions very close to 
coordinate singularities, where the source terms might diverge. 



4. Tests and applications. 

In order to demonstrate the feasibility of our procedure we have considered an exhaustive 
sample of standard discontinuous initial value problems for which the exact solution is 
known, as well as some numerical applications involving strong gravitational fields. The 



set of SRRS used in the computations are the exact one (Marti and Miiller 1994) for ID 
problems, as well as SR extensions of the linearized solvers described in Harten, Lax and 
van Leer ( |1983[ ), Roe ( |1981| ), and Donat and Marquina Ql996| ). 

To summarize, we have successfully redone all the experiments shown in Banyuls et 



al. (1997), that includes relativistic shock-tube tests for non-diagonal metrics, as well 
as a number of simulations of relativistic wind accretion onto a Schwarzschild black 
hole. In Figure 1, we show the results from a simulation of spherical accretion of an 
ideal gas onto a Schwarzschild black hole. The analytical solution derived by Michel 



(1972) is represented by the solid line and the diamonds represent the numerical solution 
obtained using the exact SRRS, after the stationary state has been reached. In Figure 
2, we display results from two dimensional simulations of non-spherical accretion onto a 
moving black hole, corresponding to one of the models recently studied in Font & Ibanez 



( |l998| ) [model MC2 in their table I]. The figure displays isocontours of the rest mass 
density in logarithmic scale. The upper-left panel displays the results obtained from the 



code described in Banyuls et al. (1997), the upper-right panel shows the results obtained 
with the new approach using the same solver (Roe- like). Results obtained with the new 
approach using HLLE and Marquina's solvers are shown in the lower-left and lower-right 
panels, respectively. 

The main conclusion emerging from the comparison is that our new approach gen- 
erates remarkably similar results for the four different SRRS, the tiny differences being 
due to the intrinsic properties of each solver, e.g., Roe's solver is the least diffusive and 
therefore more oscillatory. The results following the method presented in Banyuls et al. 
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l997[ ) are indistinguishable from the ones obtained using the special relativistic Roe- 
like solver in this new approach. It can be shown analytically that both algorithms are 
equivalent. 

5. Conclusions and outlook 

We have developed a general procedure to use SRRS in multidimensional general rela- 
tivistic hydrodynamics that allows us to take advantage of the increasing number of SRRS 
developed recently, overcoming partial approaches derived in previous papers, which were 
restricted to linearized Riemann solvers. Since the change of coordinates we propose is 
linear and only involves a few arithmetical operations, the additional computational cost 
of the approach is completely negligible. 

The procedure has a large potentiality and could be applied to other systems of con- 
servation laws, as the equations of general relativistic magneto-hydrodynamics, providing 
a very useful numerical tool to solve them using the corresponding Riemann solvers devel- 
oped for the special relativistic case. The feasibility of the approach has been extensively 
tested and its numerical performance is, at least, as good as other schemes developed 
in previous papers, having the additional advantage of being very well suited to include 
future work and improvements that might be done in the field of SR Riemann solvers. 
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Fig.l. Spherical accretion of an ideal gas: profiles of density (p), internal energy (e) and 
velocity (v/c) as a function of the radial coordinate, after the steady state has been reached. 
The critical point (r c ), the critical density at the critical point (p c ) and the adiabatic exponent 
of the equation of state (7) have been taken r c = 200M, p c = 7 x 10~ 4 and 7 = 4/3. The solid 
line corresponds to the analytic solution and the diamonds to the numerical one obtained using 
the exact SRRS. 

Fig. 2. Non-spherical hydrodynamic accretion onto a Schwarzshild black hole. The initial 
model is characterized by an asymptotic velocity of 0.5c and Mach number 5. The adiabatic 
exponent of the fluid is 5/3. The simulation employs a grid of 120 x 40 zones in the radial and 
polar directions respectively. The radial domain extends from 2.1M to 38M. The polar domain 
extends from to it. The flow moves from left to right. The different pannels show isocontours 
of the logarithm of the density normalized to its asymptotic value. Starting from the upper-left 
pannel and in a clockwise sense we show results for the Roe solver (as used in Banyuls et al 
(1997)), its SRRS version, Marquina's solver and HLLE. The isocontours span the same interval 
regardless of the solver used. This range goes from to 1.15. The maximum values are always 
found at the rear part of the hole where the matter piles up. The different solvers agree on this 
value up to three significant figures. 
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